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Abstract 

We propose an unbiased Monte-Carlo estimator for - ,W„)], where X 

is a diffusion process defined by a multi-dimensional stochastic differential equation 
(SDE). The main idea is to start instead from a well-chosen simulatable SDE whose 
coefficients are updated at independent exponential times. Such a simulatable process 
can be viewed as a regime-switching SDE, or as a branching diffusion process with one 
single living particle at all times. In order to compensate for the change of the coeffi¬ 
cients of the SDE, our main representation result relies on the automatic differentiation 
technique induced by Bismu-Elworthy-Li formula from Malliavin calculus, as exploited 
by Eournie et al. [14] for the simulation of the Greeks in financial applications. In 
particular, this algorithm can be considered as a variation of the (infinite variance) 
estimator obtained in Bally and Kohatsu-Higa [3, Section 6.1] as an application of the 
parametrix method. 

MSC2010. Primary 65C05, 60J60; secondary 60J85, 35K10. 
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1 Introduction 

Let d > 1, T > 0 and W be a d-dimensional Brownian motion, /x ; [0,T] x —)• 

and cr : [0, T] X — )• be the drift and diffusion coefficients, where denotes 

the collection of all d x d dimensional matrices. Under standard assumptions on these 
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coefficients, we consider the process X defined as the unique strong solution of the 
multi-dimensional SDE, 

Xo = xo, and dXt = fi[t,Xt) dt -|- a{t,Xt) dWt, (1.1) 

Our main focus in this paper is on the Monte-Carlo approximation of the expectation 

Fo := (1.2) 

for some function g : —>• M and discrete time grid Q < t\ < ■ ■ ■ < tn = T. 

When n = 1, the analytic formulation of the problem is obtained by the well-known 
representation Vo = 'w(0,Xo), where u is the solution of the linear PDE 

diU + b{t,x) ■ Du + ^^T\ao'"'{t,x)D‘^v\ ut = g, (1.3) 

In practice, the Monte-Carlo method consists in simulating N independent copies of 
a discrete-time approximation of X, and then estimating Vq by the empirical mean 
value of the simulations. The corresponding error analysis consists of a statistical error 
induced by the central limit theorem, and a discretization error which induces a biased 
Monte Carlo approximation. Under some smoothness conditions, Talay and Tubaro 
[25] proved that the discretization error for the Euler scheme is controlled by a rate 
CAt, where At denotes the time step discretization. Since then, many works focused 
on the analysis of the discretization error under various discretization techniques, see 
e.g. Kloeden and Platen [23], and Graham and Talay [17] for an overview. However, 
the statistical error estimate N ~2 is lost in all cases, as its combination with the 
discretization error leads to an overall error estimate of the order for some 

e > 0. 

In the context of one-dimensional homogeneous SDEs with constant volatility coef¬ 
ficient, Beskos and Roberts [6] developed an exact simulation technique for X by using 
the Girsanov change measure together with a rejection algorithm, see also Beskos, Pa- 
paspiliopoulos and Roberts [7], Jourdain and Sbai [22], etc... This technique also 
applies to more general SDEs by using of the so-called Lamperti transformation which 
reduces the approximation problem to the unit diffusion case. We also refer to the 
subsequent active literature of exact simulation of an L°°—approximation of X, see. 
e.g Blanchet, Chen and Dong [5]. 

An alternative approximation method for Vq was induced by the multilevel Monte- 
Carlo (MLMC) algorithm introduced by Giles [15], which generalizes the statistical 
Romberg method of Kebaier [18]. One of the main advantages of the MLMC algo¬ 
rithm is to control the global error (sum of discretization error and statistical error) 
with a much better rate w.r.t. the computation complexity. We refer to Giles and 
Szpruch [16], Alaya and Kebaier [1], Rhee and Glynn [24] for further developments. 
In particular, Rhee and Glynn [24] obtained an unbiased simulation method for SDEs 
by using a random level in the estimator. 

More recently. Bally and Kohatsu-Higa [3] provided a probabilistic interpretation 
of the parametrix method for PDEs. In particular, when n = 1, they obtained a 
representation formula for Vq in form 

E[5(At)>Vt], (1.4) 
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where X is defined by a Euler scheme of X on a random discrete-time grid (the 
time step follows an independent exponential distribution), and Wt is a corrective 
weight function depending on X. The above representation is formally similar to 
the stochastic finite element method proposed by Bompis and Gobet [8], where one 
replaces X by its Euler scheme solution in (1.2) and then corrects partially the error by 
some well-chosen weight functions. Notice that in the above representation of Bally and 
Kohatsu-Higa [3], the process X can be exactly simulated and hence it may provide an 
unbiased estimator for Vq. Nevertheless, the obtained weight function Wt is integrable 
but has infinite variance, and hence the corresponding Monte-Carlo estimator looses 
the standard central limit error estimate. 

In this paper, we provide a representation of Vq in the spirit of (1.4), but with an 
alternative weight function for the representation. Our results follow from completely 
different arguments. More importantly, our unbiased approximation of Vq has finite 
variance, and applies for a large class of SDEs. 

Our main idea is to consider the Euler scheme solution X as solution to a regime¬ 
switching SDE with some well-chosen coefficients. In order to compensate for the 
change of the coefficients of the SDE, we introduce some weight functions obtained by 
the automatic differentiation technique induced by Bismut-Elworthy-Li formula from 
Malliavin calculus, as exploited by Eournie et al. [14] for the simulation of the Greeks 
in financial applications. 

The technique introduced in the present paper is inspired by the numerical algo¬ 
rithm introduced in [19, 21], for semilinear PDEs of the form 

dtu + ^Au + Fo{t,x,u) =0, UT = g, 

for some nonlinearity Eq. The main idea in [19, 21] is to use an approximation by a 
branching diffusion representation induced by approximating the nonlinearity Fq by a 
polynomial in u. Namely, given the nature of the linear operator, the representation 
is obtained by means of a Brownian motion with branching driven by the polynomial 
approximation of Fq. 

Loosely speaking, the method developed in the present paper follows by reading 
the PDE part of (1.3) in the following equivalent form: 

dtu +-Au + Fi{t,x, Du, D'^u) = 0, 

where 

Fi{t,x,z,j) := b{t,x) ■ z + ^Tr[{aa'^{t,x) - I)j]. 

However, in contrast with the nonlinearity Fq, the above function Fi involves the 
gradient and the Hessian of the solution u. Consequently the last PDE can not be 
handled by the existing literature on branching diffusion representation of PDEs. The 
automatic differentiation technique introduced in the present paper is an important 
new idea which allows to convert Du and D^u in Fi into u. Since no powers of u are 
involved in the equation, this leads to a representation by means of a Brownian motion 
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with exactly one descendent with two different possible types revealed by the weight 
function corresponding to the order of differentiation. 

We believe that the automatic differentiation trick introduced here has very im¬ 
portant consequences, beyond the particular application of the present paper. Indeed, 
in our paper [20], we provide a significant extension of the branching diffusion repre¬ 
sentation to a general class of semilinear PDEs. 

The paper is organized as follows. In Section 2, we consider the SDE with constant 
diffusion coefficient, and propose an unbiased estimator for Vq for both Markovian 
case and path-dependent case. Then in Section 3, we consider the SDE with general 
diffusion coefficient function, and obtain a similar representation formula for Vq, which 
is integrable but of infinite variance. Next, in Section 4, we provide some numerical 
examples. Finally, we complete some technical proofs in Section 5. In particular, 
an easy example is studied in Section 5.1 to illustrate the main idea of the technical 
proofs. 


2 Unbiased simulation of SDE with constant 
diffusion coefficient 

In this section, we will restrict to the constant diffusion coefficient case, and propose 
an unbiased estimator for Vq having finite variance. 


2.1 The Markovian case 


Let us start by the Markovian case, where the diffusion process X is defined by 


Xo = xo, dXt = fi{t,Xt)dt + ao dWt, (2.1) 

for some matrix co S M'^, and our objective is to compute 


Do = ngiXr)]. 


( 2 . 2 ) 


for some function g : —)• M. We impose the following conditions on /i and (Tq: 


Assumption 2.1. The diffusion coefficient ao is non-degenerate, the drift function 
g,{t, x) is bounded continuous in {t, x), uniformly ^-Holder in t and uniformly Lipschitz 
in X, i.e. for some constant L > 0, 


g.{t,x) - fi{s,y) 


< L(^y/\t - s| -h |x - y|^, V(s,x), (t,y) e [0,r] X M'^. 


(2.3) 


2.1.1 The unbiased simulation algorithm 

To introduce our unbiased simulation algorithm, let us first introduce a random discrete 
time grid. Let /3 > 0 be a fixed positive constant, (rj)j>o be a sequence of i.i.d. T(/3)- 
exponential random variables. We define 
k 

Tfc := A: > 0, and W := max {A: : < t|. (2.4) 

i=l 
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Then {Nt)o<t<T is a Poisson process with intensity j3 and arrival times (Tfc)fc>o- We 
denote also Tq := 0 and AT^+i := T^+i — Tk- 

Let VP be a d-dimensional Brownian motion independent of (rj)i>o, we introduce 


:= Wt^ — Wt^_^, k > 0. 


and a process X as the Euler scheme of X on the random discrete grid (rfc)fc> 0 ) i-e- 
Xq = xq and 

^Tk+i ■= ^Tk + k'{Tk,XT^)ATk+i + croAVVTi,_^j, A: = 0,1,-- - ,Nt. (2.5) 
Then our estimator is given by 


— ^/3T 


if} := e 


g (^Xt^ - g ) l{ArT>o} 


Nt 

/3-iVT 

k=l 


ki 


( 2 . 6 ) 


with 


1 (/i(Tfc, AtJ - /i(Tfc_i, Ar,_,)) • (4)-^A1 Pt, 


W: := 


fe + 1 


ATfc+i 

Theorem 2.2. Suppose that Assumption 2.1 holds true, and g is Lipsehitz. Then 

]E[(^^] < oo and Vb = 1 E[ ^]. 

v 2 n 


(2.7) 


Proof, (i) We first show that E[('?/’) J < cc- For simplicity, we denote AA^ := 
Xxf. — ^Tfc_i for A; > 0. Let Lg be the Lipsehitz constant of the function g, and set 
Lq := |(o'oO'([) > 0 by the non-degeneracy of do. Then using Assumption 2.1, it 

follows by direct computation that 

Nj' '' 


e-^^^\<Lg(\g{xo)\+AT, + \AXTA) H 


L(v/M^+|AAt,+J) 


k=l 


PATk+i 


(4)-^AWt, 


Then denoting := E[ • | AT^+i], we have 


v^MuT + |aa 


T’fc-l-ll / T\-l 


ATi 


k-\-l 


(d(J )-^AVTt, 


fc +1 


2 i 


< E 


fe+i 


{l + \f,\ooVT+\aoZ\f\{a^)-^Z'^ 


where \g\oo ■= ^jYli=i \hi\o^ Ik'ilo •= sup^ 3 , |/rj(t, x)|, and A is a standard centered 
normal distribution in This provides 




\/AFfc+i -|- |AAj'j.+i| 


ATi 


k+l 


[alr^AWT,^^ 


< 2 (l + |/i|ooVr)2 E[|(d([)-'A|'] + 2 E[|doA|'|(d([)-iA|'] 
= 2(1 -V |/i|ooA/T)^ Tr((crocr([)"^) -V 2{3d + d{d-l)) =: 7 . 


We therefore get the following upper bound: 

E['02] < , where C := L 2 e[(| 5 ( 0 )| + ATi + |AAi|)^]. (2.8) 

(ii) The equality Vq = E[^] will be proved in Section 5, with illustration of the main 
idea in Section 5.1. □ 
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2.1.2 On the choice of /5 


Notice that the random variable if) in (2-6) can be exactly simnlated from a seqnence 
of Gaussian AA(0,1) and exponential £{I3) random variables. Then the integrability 
and representation results in Theorem 2.2 induce an unbiased simulation Monte-Carlo 
method to approximate Vb, with error induced by the standard central limit theorem. 

We next observe that the constant /3 > 0 may be chosen so as to minimize the 
approximation error relative to the computational effort: 

• By the central limit theorem, the error induced by the Monte Carlo estimator 
based on the representation ip is characterized by the variance of ip. For tractabil- 
ity reasons, we shall instead replace it by the bound (2.8). 

• The computation effort is proportional to the number Nt of arrivals of the Poisson 
process before the maturity T, and is thus given by C'K[Nt] = CfiT. 

In view of this, we shall choose /3 by minimizing the ratio of the variance bound (2.8) to 
the mean computational effort. This minimization problem is obviously independent 
of the constants C, C, and reduces to: 

mm/(/3), where /(/3) := ^ exp (^r(/3 -h ^))- 

Direct computation shows that the equation /'(/3) = 0 has a unique solution on (0, oo) 
given by 

/?* := ^ 7 ^ 2 + 2 - 2/4 + 'L, 

As lim^^o/(/3) = hm^^oo/(/3) = oo, this shows that (5* is the minimizer of the 
above defined criterion, and will be taken as our “best sub-optimal” choice of /3 for 
the unbiased estimator ip. 

2.2 The path-dependent case 

In this part, we would like to provide an extension of the above estimator ip in (2.6) 
to the path-dependent case. Let n > 0, 0 = to < ti < ■ ■ ■ < tn = T, ao ^ be a 
non-degenerate matrix, and /i : [0,T] x —>• be a continuous function, Lipschitz 

in the space variable. Let X be the unique solution of SDE, with initial condition 
Aq = xo, 


dXt — Xti/\tr " ^ ^t„At) dt + aodWt; (2.9) 

and the objective is to compute the value, 

Vh := E[g{Xt„--- (2.10) 

for some Lipschitz function g : M. 


6 




Remark 2.3. It is clear that the value Vq defined above can he characterized by a 
parabolic PDE system. Namely, for every k = 1, - ■ ■ ,n and (xi, • • • , Xk-i) S 
we define 

hk{t,x) := fi{t,xi,--- ,Xk-i,x,--- ,x), V(t,x) e [4_i,4] X (2.11) 

Suppose that {uk)k=i,- ,n is a family of functions such that Uk is defined on [tk-i,tk] x 
^dxk and x Uk{t,xi, ■ ■ ■ ,Xk-i,x) is a solution (at least in the viscosity sense) of 

dtUk + ^CToo-Q : + hk ■ Duk = 0, ( 2 . 12 ) 

with terminal conditions 

^k {ik ) ) ■ ■ ■ ) ®fc) ^fc +1 {I'k ) j ■ ■ ■ ) j ^fc)) k I)'"' 1 ) 

and Un{tn, xi, - ■ ■ , Xn) = g{xi, • • • , Xn). Then we have Vq = iti(0, xq). 


2.2.1 The algorithm 

The unbiased simulation algorithm of Vq can be obtained by an iteration of the estima¬ 
tor (2.6) on each time interval [tk, tfc+i]. One should just be careful on the integrability 
issue. Let us first introduce the algorithm. 

Recall that IT be a standard d-dimensional Brownian motion, (rj)j>o is a sequence 
of i.i.d. £l(/3)-exponential random variables independent of W. Then N = {Ns)Q<s<t 
and {Ti)i^Q are defined in (2.4). Define further for every k = 1, • • • ,n, := Nt^, — 

Nt^_^ the number of jump arrivals on [tk-i,tk), and Tq := tk-i and Tj := 

ATf := Tf - T/_i, := Wfk, NWf := Wf - Vj = 1, • • • , iV^ + 1. 

Example 2.4. We give below an example for the case n = 2. In the following example, 
the number of jump arrivals on [ 0 ,ti) is = 2, that on [ti,t 2 ) is = 1, and total 
number of jump arrivals is Nt = 3. 

For k = 1, we have Tq = 0, T) = Ti, = T 2 and T^ = h; Wq = 0, W( = VTti; 
W) = Wt 2 and Wq = Wtj^. For k = 2, we have Tq = ti, Tf = Tq, T 2 = t 2 , and 
= Wt^, Wf = Wt 3 and = Wt^- 


Ti 


T2 


ti 


Tq 


t2 


We next introduce a process , Vj = 0,1, • • • , Nk + 1, for each k = 1, ■ ■ ■ ,n 

and initial condition x = {xq, xj, - ■ ■ , Xk-i) S by Xq’^ := Xk-i and 

rk,'X.\ 


^ + croAWjVi. 


Similarly, for every j = 1, • • • , Nk, we define a automatic differentiation weight, with 
Uk defined by (2.11), 

Wi)) ■ 


W 

3 


^^j+l 
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We now introduce the algorithm for the path-dependent case, in a recursive way. 
First, for x = {xo,xi,--- ,Xn) G set i>n+i ■= '' ,Xn)- Next, for 

k = 1, - ■ ■ ,n, denote 

:= ,Xk-i, and := (xq, xi, • • • ,Xfe_i, 


Then given V'fc+u we dehne 

■= e' 






k+l V'fc+i -*-{Ar''>0} 


/ 3 "^ ( 2 . 13 ) 

j=i 


We finally obtain the numerical algorithm of the path-dependent case: 

^ := 


( 2 . 14 ) 


2.2.2 The integrability and representation resnlt 

We notice that the algorithm in the path-dependent case is nothing else than an 
iterative algorithm of the Markovian case, as suggested by the PDFs (2.12) in Remark 
(2.3). When the random variable 'ip in (2.14) is integrable, it is not surprising to obtain 
the representation Vq = E[t/’] as a consequence of Theorem 2.2. However, because of 
the renormalization term (i.e. ^{7V'=>o}) ™ (^•^2))! variance analysis 

becomes less obvious. We provide here a sufficient condition to ensure that ip has finite 
variance. 

Theorem 2.5. Suppose that p, : [0,T] x —>• and g : —)• M are differen¬ 

tiable up to the order n, with bounded derivatives. Then 

^[("0)^] < oo and Vq := ^[-0]. 

We will prove the integrability result here, and leave the proof of the representation 
result Vq := E[V’] in Section 5. As preparation, let us first provide two technical 
lemmas. Let tt = (0 = sq < si < • • • < Sm = T) be an arbitrary partition of 
the interval [0, T], p : [0,T] x —>■ a valued function. We define by 

:= X and 

+ p{sk,Xl’^)Ask+i + Ws,^,-Ws,. (2.15) 

Further, let —>• M be a smooth function, i > 0 and i = (ii, - ■ ■ ,i£) G {!,■ ■ ■ , dY, 

we denote di^ifix) := di.^...^.^ip{x). 

Lemma 2.6. Suppose that x >->• p{t, x) is differentiable up to order n with uniformly 
bounded derivatives, and X'"'^ is defined by (2.15) with initial condition Xq’^ = x. 
Then x i—)■ X^^ is differentiable up to order n and there is a constant C independent 
of the partition n such that 

max max max < C. 

l<e<n 0<fc<m I « I 


Proof. For simplicity, we consider the one dimensional d = \ case, while the multi¬ 
dimensional can be deduced by almost the same arguments. First, let = 1, we 
have 

+ d,fl{sk,X'^^’^)d^Xl'^Ask+i, 

which implies that 

= U']tl{l + d,fl{sk,Xl’^)Ask+i). 


Since dxfl{t, x) is uniformly bounded, it follows that dxX^' is bounded by some con¬ 
stant Cl independent of 1 < A: < m and the partition tt. By induction, it is easy to 
deduce that for £ = 2, • • • , n. 




where Pi is a Polynomial on d^^ip,{sk,X^’^) and for i = 1, • • • ,i — 1, which is 

uniformly bounded by some constant independent of A: = 1, • • • ,m and the partition 
TT. Hence d^gX^’^ is also bounded by some constant Ce independent of A: = 1, • • • , m 
and the partition tt. □ 

Lemma 2.7. Let {'ip^)i<k<n+i defined by (2.13). Then for every k = 2, - ■ ■ ,n -|- 1, 
and every x = (xo,xi,--- ,Xk-i) G the map Xk-i derivatives up to 

order k — 1 and 


max 

i<e<k-i 




< cYliN^ + iy-\ 

j=k 


(2.16) 


Proof. We will prove it by induction. First, let k = n-|-l, then fin+i ■— xi, - ■ ■ , x„) 
and hence < C for some constant C and for every £ = 1, • • • , n. 

Next, suppose that (2.16) holds true for we know from (2.13) that 




xfe.x 

'fc-i-i 


'T'vfciX.O _ 

- ^f+i 1 


N 


{iV'“>0} 




X 


k,x. ' 


j- 1 ^ 


K) 




j+i- 


Then using the estimation in Lemma 2.6, we see that (2.16) is also true for and 
we hence conclude the proof. □ 

Proof of Theorem 2.5 (i). By Lemma 2.7, we know that x i—>■ is differentiable 

and in particular uniformly Lipschitz with coefficient bounded by 2Cn”^2(-^'^ + 1)'^~^- 
Then the definition of falls into the Markovian case n = 1, but with terminal 
condition x i—>■ ■ Notice that < Nt has a Poisson distribution: ¥{Nt = m) = 

g-/3r (/3T) ^ follows that, for some constant C > 0, 


E[|V’r“|^] < E 


JJ(W + 1 ) 

J=2 


2(i-i) 


< E 


4^2^7Vr^^^^l)n(n-l) 


< oo, 


which implies that has hnite variance. 


□ 
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3 Unbiased simulation of general SDEs 

Let us now consider the SDE (1.1) with general diffusion coefficient function, i.e. with 
drift and diffusion coefficients fi : [0, T] x —)• and a : [0, T] x —>• 

Xq = xq, and dXt = fj,(t,Xt) dt + a(t,Xt) dWt- 


Our objective of study in this section is 

Vo = E[ 5 ((X 7 ’)] , for some function <7 : —)• M. 


We will provide a representation result of Vq in the same spirit of that in Section 2. 


Remark 3.1 (Lamperti’s transformation). We also notice that in some cases, the 
above SDE (1.1) may be reduced to the constant diffusion coefficient case (2.1), by the 
so-called the Lamperti transformation. 

(i) When d = 1 and a{t,x) > 0, let us define a function /i : [0, T] x M —>■ M 

rx j 

h{t,x) := / dy. 

Jo cr{t,y) 


Notice that for fixed t € [0, T], x >->■ h{t,x) is strictly increasing, we denote h~^(t,-) 
its inverse function. Then by ltd’s formula, it is easy to obtain that Yt := h{t,Xt) 
satisfies the SDE 


dYt 


dth(t,h ^{t,Yt)) + 


fi{t,h ^{t,Yt)) 

a{t,h-^{t,Yt)) 



dt + dWt, 


whose diffusion coefficient is a constant as in SDE (2.1). 

(ii) When d > 1, a is non-degenerate and satisfies some further compatibility conditions, 
one can also obtain a similar transformation to reduce SDE (5.8) to the constant 
diffusion coefficient case. 


3.1 An estimator of infinite variance for general SDEs 

Let us impose the following conditions on coefficient functions p. and a. 

Assumption 3.2. The function {fi,cr) : [0,T] x x and a := ^crcr"'' : 

[0, T] xM'^ —)• are uniformly bounded, and are uniformly Holder in the time variable, 
uniformly Lipschitz in the space variable, i.e. for some constant L, 

\{fj,,a,a){t,x) - {iJ.,a,a){s,y)\ < L{y/\t - sj + |x - y\), (3.1) 

for all {t,x),{s,y) G [0, T] x and a{t,x) is non-degenerate such that, for some 
constant Eq > 0, 

a{t,x) := ^aa^{t,x) > Sold, V(L a:) G [0, T] x M'^. 
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Recall that {Tk)k>o are defined by (2.4) with a sequence of i.i.d. £l(/3)-exponential 
random variables, and bb is a Brownian motion; the increment of the Brownian motion 
are defined by := bbt^. — bbt^_^, and AT^ := — T^-i- As in (2.5), we introduce 

X as solution of the Euler scheme on discrete grid by Xq := xq and 

^Tk+i ■= XTk+fJ-{Tk,XTi^)ATk+i + a{Tk,XT^^)AWT^^-^ /c = 0, • • • , At- (3.2) 


We then introduce a representation formula by 

where, for each A; = 1, • • • , Nt, 


Nt 


k=l 


r ^ ^ ^ M {^'^iTk,XT,))-'AWT, 

■= [KTk,XT^) - fi{Tk-i,XT^_^)\ ■ - ATk+i - 


fc + 1 


and 


Wl := [a{Tk,XT,)-a{Tk-i,XT,_,)] 


{a^{Tk,XT,)) 


_i AWT,+,AbbT^^ - ATk+,Id 


(3.3) 




a{Tk,XT,)-^ .(3.4) 


Theorem 3.3. Suppose that Assumption 3.2 holds true, and g is Lipsehitz. Then 


E[|^|] < oo and Eo = lE[^/^]. 

Proof, (i) Consider the random vectors := and := --, for 

all k = 1, • • • ,Nt + 1, which are independent of ATj. conditional on {AT^ > 0} = 
{Nt > k — 1}, and which have finite second order moment. Notice that g{t,x) and 
a{t, x) are uniformly bounded, and 1/2—Holder-continuous in t and Lipsehitz in x, and 
a is uniformly bounded from below above zero. Then, for each fe = 1, • • • , Nt, 


\wl\<c[^k + \Xn-XT,_,\) 


ATfc+i 


< C(l -I- v^ATfc -b 



where the constant C > 0 may vary from term by term but is uniformly bounded for 
all k. Similarly, one obtains that 


\wl\ < C(l + v^+ldl)ld+i 

As ATk < T, it follows that 


ATk 


ATk+l 


\wl\ + \wl\ < C{1 + le^l) {la+A + Id+ll) 

for some constant C > 0 independent of k. In addition, we have by the Lipsehitz 
condition on g that 


E 


\9{Xt) - 9{Xtj.^ )\ 


AT, 


A^j’+l 


< C\J AT/Vy+i. 
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Then, it follows from the expression of in (3.3) that 

N'j' 


E[l<i|] < CE[n 

k=l 
Nt 

^[n 


c 


Nt 


< CE 


C 




fc +1 




+ CE 


+ CE 


\g{XT)\ l{AfT=o} 

\g{xo + /r(0,xo)r + cr(0 , xo)Wt)| 


for some constant C > 0, where we have also nsed the independence of the ^^’s and 
their the boundedness of their second order moments. The integrability of ip is now a 
direct consequence of Lemma A.2. 

(ii) The proof of the equality Vq = E['i/)] will be completed in Section 5. □ 

To conclude, we notice that the variable pj is of order in general 

cases, and the latter is integrable but of infinite variance. Therefore, ip is not a good 
estimator for Monte-Carlo method. Nevertheless, it should still have some theoretical 
value as an alternative representation formula obtained by Bally and Kohatsu-Higa [3, 
Section 6.1]. 


3.2 An estimator for one-dimensional driftless SDE 


To overcome the problem of variance explosion of the estimator (3.3), we will consider 
the higher order approximation X of X, and obtain an estimator of finite variance for 
the one dimensional (d = 1) driftless SDE of form 

Xo = xo, dXt = a{t,Xt) dWt, (3.5) 


Our objective is to compute 


Vo := E[g'(X'r)], for some function s': M—)• M. 


Recall {Tk)k>o has been introduced in (2.4) from a sequence of i.i.d. exponential 
random variables, independent of the Brownian motion W. We next define X by 
-^0 = Xq, 

dXt = (a{Tk,XT,) + dMTk,XT,){Xt - XtJ) dWt, on [T^,T^+i], (3.6) 

for A: = 0,1, • • • , Nt- By denoting 

4:=aiTk,XT,)-dMTk,XT,)XT, and 4 ■= d^Tk, Xt,), (3.7) 


then the above linear SDE (3.6) has an explicit solution which is given by 
XT,^,=XT,+a{n,XT,)XWT,^„ if c^ = 0, 

and 


(3.8) 


Xt, 


k + 1 


Cl C'l / 

= - ^ ^ exp - 


r-k\2 


-ATfe+i + c^AWt, 


A:-1-1 


+ Xt, exp ( - Ml^ATfc+i + c^AWt, 


k + 1 / ’ 


if 4 + 0. (3.9) 
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We then define ijj by 


^ := 


giXr) - 5(^riVj,)l{Arr>0} 


/3 


N'T 

—Nt T T ttt2 




(3.10) 


k=l 


where the automatic differentiation weight is given by (see Lemma 5.9 below) 


yy2 _ a{Tk,XT^)-ak 


2a{Tk,XT^. 


- dMTk,XTj 


, AWl^^-ATk+,^ 


ATi 


+ 


k-\-l 




J, (3.11) 


with a(-) := ak := and dk ■= cT(Tfc_i,J + dx(T{Tk-i,XTk-i){XTu: - 

Xn-i)- 

Similarly to the discussion at the end of Section 3.1 (see also Remark 5.8 below), 
the variable ip in (3.10) is integrable but of infinite variance in general. To make the 
variance finite, we introduce an alternative estimator using an antithetic variable. Let 
X^ be an antithetic variable of Xt defined by 


:= - a{TM^,XT^)AWT^, if = 0 , 


and 


y- — _ 2l _ 

- Nt 


+ 


_ 

Nt 


exp 


„^t\2 


-ATy^+i - C^^AWt, 


NrTi-\-l 


+ Xt^^ exp ( - - c^^AWt, 


Nrp+1 J 5 


if / 0 . 


Denote Wk ■= for A; = 1, • • • , Nt — 1 and 


W 


Nt 


(x{Tnt ; XJTj, ) Q^Nt 


2Qi(TNt 7 XjTrp 

We then introduce 

"0 + V’ 


^x^^^Nt 7 XJTj- 


AWt^ 


AWl 


AT; 


+ 


TJV 7.+1 


- AT, 








N’j' 


'll; := 


with g{X^) - g{XT^^)l{y^^o} H ^*7 • (3-12) 


k=l 


Notice that the Brownian motion is symmetric, thus ip has exactly the same distri¬ 
bution as 0 , and it serves as an antithetic variable. 

Assumption 3.4. The diffusion coefficient cj(-) satisfies a{t, x) > e > 0 for all (f, x) G 
[0,T] X M, a{t,x) is bounded and Lipschitz in {t,x), dxcr{t,x) is bounded continuous 
in {t,x) and uniformly Lipschitz in x. Further, the terminal condition function gf) E 

cim- 

Theorem 3.5. Suppose that Assumption 3-4 holds true. Then 

E[|^|] + E[|^|^] < 00 ; and Rq = E[^ = E[^. (3.13) 

We will complete the proof in Section 5.4. 
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Remark 3.6. (i) As if: has finite variance, we may use the representation of The¬ 
orem 3.5 to built an unbiased Monte-Carlo estimator of Vq. However, given the as¬ 
sumed regularity conditions, and the restriction to the one-dimensional setting, such 
a Monte-Carlo approximation is not competitive with the corresponding PDE based 
approximation methods. However, we believe that the present methodology is open to 
potential improvements, and we hope to improve our results in some future work so as 
to address the higher dimensions. 

(ii) For a general SDE with drift function and/or d>l, we can also consider a similar 
choice of which leads to fi{t,x) = ci + C 2 X and a{t,x) = 03 + c^x and a linear 

SDE 


dXt = {dC2Xt) dt + {c3-\-C4Xt)dWt, (3-14) 

where ci £ C 2 , C 3 € and C 4 is linear operator from to M'^. However, to the 

best of our knowledge, the exact simulation of linear SDE (3.14) in high dimensional 
ease, as well as the associated automatic differentiation (Malliavin) weight as in (3.11) 
(see also Lemma 5.9 below), is still an open question. 


4 Numerical examples 

Notice that our estimator if given by (2.6) (resp. if given by (2.13) and (2.14)) is an 
unbiased estimator for Vq in (2.2) (resp. Vq in (2.10)). Then the error analysis of the 
Monte-Carlo approximation reduces to the statistical error. Hence the computation 
cost to achieve the accuracy 0 (e) for the approximation of Vq (resp. Vq) is of order 
0 (e“^), thus avoiding of the dependence on the discretization error. 

By combining different level of simulations, the MultiLevel Monte Carlo (MLMC) 
method proposed by Giles [15] achieves a computation cost of order 0(e“^(loge)^) or 
0(e“^) depending on the strong discretization error rate. In particular, by considering 
a randomization of the level, Rhee and Glynn [24] obtained an unbiased estimator. 
In the following, we provide some numerical results and comparisons between our 
unbiased simulation method with the Euler based MLMC method proposed by [15]. 

4.1 Two one-dimensional SDEs 

Let VL be a one-dimensional standard Brownian motion, we consider the SDE given 
by 

5o = l, dSt = O.liy'M ASt-l)Stdt + ^StdWt, 

where M is a large constant introduced in order to guarantee the Lipschitz property 
of the drift coefficient (in our numerical implementation, we have observed that the 
value of M is not relevant for large M, and that the numerical finding are not changed 
by taking M = 00 ; this hints that our results may be extended beyond the case of 
Lipschitz coefficients). Applying Lemperti’s transformation Xt := log(5i), we reduce 
the above SDE to the constant diffusion coefficient case, in form of (2.1), 

Ao = 0, dXt = (o.1{VM a - l) - l/s) dt + ^dWt. (4.1) 
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Mean value 

Statistical error 

Computation time 

US {N = 10®) 

0.204864 

0.00140709 

0.016814 

MLMC 

0.204993 

0.000949166 

0.032017 

US {N = 10®) 

0.205396 

0.000444462 

0.171835 

MLMC 

0.205602 

0.000308634 

0.234526 

US {N = 10^) 

0.20552 

0.000142554 

1.63013 

MLMC 

0.205648 

0.0001 

1.96197 

US {N = 10®) 

0.205641 

4.52282e-05 

16.2189 

MLMC 

0.205638 

3.18855e-05 

18.3833 


Table 1: Numerical results for Vq in (4.2) (case d = 1), US denotes our unbiased simulation 
algorithm with /? = 0.1, the computation times are expressed in second. 


We implement our unbiased simulation method for the two following expectations; 


Uo := E[(5 t- 1^)+], and Uq 




(4.2) 


where we choose K = 1, T = 1, n = 10 and ■= ^T. Notice that the path- 
dependent example does not satisfy the differentiability sufficient condition in Theorem 
2.5. However, our numerical findings do not show any numerical difficulty in the 
present setting. 

Using different numbers N of simulations, we obtain the standard deviation as 
(statistical) error of our estimator. Next, using the errors obtained by our unbiased 
simulation method, we implement the MLMC algorithm in Section 5 of Giles [15], and 
we compare the computation time (in second) of the two methods. More precisely, 

the statistical error of the unbiased simulation method is given by ^JYa.T['^p]/N, where 
Var[?/)] denotes the estimated variance of ?/>. For the implementation of MLMC, we 
choose M = 4, Nl = 10^ and use equation (10) in [15] as criteria to stop the loop in 
MLMC (see more details in Section 5 of [15] for the meaning of M and Nl). 

The numerical results are given in Tables 1 and 2. We observe that with the same 
Monte-Carlo error, both methods have very close performance. In the present partic¬ 
ular example, the computational time of our methods is slightly smaller. However, the 
conclusion may change depending on the nature of the example. Let us consider the 
problem 


Uo := E[sin(Xr)], 


(4.3) 


where X is defined by SDE, for some constant po £ 

Xo = 0, dXt = IJ .0 cos{Xt)dt + ^dWt- 

We implement the MLMC algorithm and our unbiased simulation method with dif¬ 
ferent value of (3, but with a given fixed error e = 0.0002. The two methods provide 
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Mean value 

Statistical error 

Computation time 

US {N = 10®) 

0.127032 

0.000762635 

0.144998 

MLMC 

0.127053 

0.000536248 

0.323337 

US {N = 10®) 

0.126363 

0.000241231 

1.40843 

MLMC 

0.126747 

0.000169842 

1.8194 

US {N = 10^) 

0.126703 

7.6418e-05 

13.9005 

MLMC 

0.126643 

5.37691e-05 

16.7499 


Table 2: Numerical results for Vq in (4.2) (case d = 1), US denotes our unbiased simulation 
algorithm with (3 = 0.05, the computation times are expressed in second. 

very close estimation of value Vq, so we give a comparison on the computation time 
in Figure 1. We can observe that /3 in the unbiased simulation method should not be 
too big nor too small, to minimize the computation effort. When po = 0.2, the com¬ 
putation time of MLMC method is slightly longer than the US method with /3 ~ 0.05. 
However, when po = 0.5, the computation time MLMC method is always smaller than 
the US method for any choice of /3 > 0. This shows that, in the context of the present 
example, the performance of our unbiased simulation method is of the order of that of 
the multilevel Monte Carlo method. 


4.2 A multi-dimensional SDE 


We next consider a d-dimensional SDE with d = 4. Let W 
4-dimensional standard Brownian motion, and (Tq the 4x4 
matrix such that 




1 

1/2 

1/2 

1/2 

1/2 

1 

1/2 

1/2 

1/2 

1/2 

1 

1/2 

1/2 

1/2 

1/2 

1 


= (W\--- be a 

be the lower triangular 


We consider the SDE 


dXt = fj.{t,Xt)dt + (TodWt, Xq = 0, i = !,■■■ ,4, 


with drift function z = I,-- - ,4) be given by fii{t,xi,--- ,Xi) = 

0.1 1 exp(xj) + |exp(a:) — 1^ — where exp(x) := (e^^ e*'‘)/4. We then 

consider two problems: 


Uo :=IE 


i 

LV4 ^ 

2=1 


and Vb := E 


1 

A4n 


k=li=l 


(4.4) 


where we choose K = 1, T = 1, n = 10 and ■= and M is a large number 
so as to ensure that the terminal condition is Lipschitz. As in the one-dimensional 
case, we implement our unbiased simulation method using different sample sizes N. 
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Comparison of the computation time of MLMC and Unbiased Simuiation method 



01_I_I_I_I_I_I_I_I_ 

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 

beta 


Figure 1: Comparison of the computation time of MLMC method and unbiased simulation 
method for problem (4.3), with the same given error. 

Then, we use the errors, obtained from our unbiased simulation method, in the MLMC 
algorithm in Section 5 of Giles [15], and we compare the computation time (in second) 
of the two methods. 

The numerical results are given in Tables 3 and 4. We observe that both meth¬ 
ods have very similar performance, with a slightly small advantage for our method. 
However, similar to the one-dimensional case, the MLMC algorithm could be better 
in other examples. 

4.3 A one-dimensional driftless SDE 

Finally, we provide an example of a one-dimensional driftless SDE. We recall that 
under the assumed regularity in Theorem 3.5, we are not expecting our method to be 
competitive with the PDE based approximations. Instead, our objective is to study 


numerically the performance of the 

estimator (3.12). 


Let us consider the SDE 




7fo = l, 

dXt = 


(4.5) 

and we aim to compute 




Fo = 

E[(Xt 

-DJ- 

(4.6) 

We implemented the simple Euler 

scheme 

with time step At = 

1/10 and simulation 


number N = 10®, and next the unbiased simulation method (3.12) with /3 = 0.1 and 
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Mean value 

Statistical error 

Computation time 

US {N = 10^) 

0.739374 

0.00921078 

0.109151 

MLMC 

0.732707 

0.00568921 

0.136884 

US {N = 10^) 

0.735745 

0.00239613 

1.06639 

MLMC 

0.733539 

0.00176862 

1.15886 

US {N = 10^) 

0.73659 

0.000831597 

10.6957 

MLMC 

0.737087 

0.000578058 

12.171 


Table 3: Numerical results for Vg in (4.4) (case d = 4), US denotes our unbiased simulation 
algorithm with /5 = 0.5, the computation times are expressed in second. 



Mean value 

Statistical error 

Computation time 

US {N = 10®) 

0.382186 

0.00247547 

0.769847 

MLMC 

0.381071 

0.00167112 

2.07589 

US {N = 10®) 

0.382846 

0.000762393 

7.65796 

MLMC 

0.383107 

0.000535905 

10.8444 

US {N = 10^) 

0.383282 

0.000244861 

85.0265 

MLMC 

0.383653 

0.00017245 

104.223 


Table 4: Numerical results for Vq in (4.4) (case d = 4), US denotes our unbiased simulation 
algorithm with (3 = 0.05, the computation times are expressed in second. 
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Mean value 

Statistical error 

Computation time 

Euler scheme 

0.161483 

0.000196733 

0.570541 

US 

0.160362 

9.34729e-05 

0.201904 

MLMC 

0.16057 

6.61696e-05 

6.65799 


Table 5: Numerical results for Vq in (4.6) (case d = 1), US denotes the unbiased simulation 
algorithm (3.12) with (3 = 0.1, the computation times are expressed in second. Notice 
also that the unbiased algorithm (3.12) contains implicitly an antithetic variance reduction, 
which makes its statistical even error smaller than that of the Eurler scheme. 

simulation number N = 10®, and then the MLMC scheme using the statistical error 
obtained from the unbiased simulation method. The results are given in Table 5, 
and we can observe that all three methods provides very similar estimation of Vq. In 
particular, the unbiased simulation method has a significant advantage. 

5 Proofs 

5.1 A toy example 

Before completing the technical part of the proofs for 
would like to illustrate the main idea by studying a 
dimensional case with unit diffusion: 

Xq = xq, dXt = fi{t,Xt)dt + dWt- 

Let 6 E M and /3 > 0, we define a sequence of i.i.d. random variable (Tfc)fc>i of 
distribution £{/3). Then let {Tk)k>i, AT^+i := T^+i -T^ and AWk+i := - ITtj, 

be defined in and below (2.4), we introduce X by 

Xt ■= xq + bt -\- Wt- 


Theorems 2.2, 2.5 and 3.3, we 
simplified example in the one 


and then define 




Nt 

5(At) 


k=l 


{^i{Tk,XT,)-b)AWT,^, 

/3ATk+i 


(5.1) 


Proposition 5.1. Let p(-,-) and g{-) be both bounded smooth functions in C"^. Then 
for all constants f3 > 0 and f3 > 0, one has 


E[|U|] < oo and K[g{XT)] = E[U]. 


Proof. Since p and g are uniformly bounded, and for some constant C > 0, the 
conditional expectation [|AWts^+i |] — V ATfc+i, then by Lemma A.2, it is 
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obvious that ijj is integral. Then it is enough to prove that ]E[ 5 (Xr)] = E['0]. In 
preparation, let us introduce 


V’n = e' 


— „/3Ti+i 


NtAu 

n 

k=l 


{^Lin,XT,)-b)AWT, 


fc + 1 


/3ATfc+i 

(^g{XT)lNT<n + dxu) (Tn+1, XT^_^_^)lNT>n^ ■ 


for all n > 0, with the convention 0^=1 = 1- is clear that (V’n)n>o are all integrable 
by Lemma A.2. 

(i) Notice that /i and g are both smooth functions, then by Feynmann-Kac formula, 
we know E[g'(X'r)] = u(0, xq), where u G C^{[0,T] x M) is a smooth function of PDE 

dtu{t,x) + -d1^u{t,x) + fi{t,x)dxu{t,x) = 0, for all (t, x) G [0, T) x M, 

with terminal condition u{T,x) = g{x). Rewriting the above PDE in the following 
equivalent way: 

dtu{t,x) + bdxu{t,x) + -dxxU{t,3o) + [fi{t,x) — b)dxu{t,x) = 0, 
it follows from the Eeynmann-Kac formula that 

rT 


tt(0,xo) = E fi((Xr) + [ {g{t,Xt) - b)dxu(t,Xt)dt 

Jo 


=/3Ti 


= E 
= E[V'o], 


9{Xt)^{Ti>t} + -^{g{Ti,XTi) - 6)9j;u(Ti, ATi)l{ri<T} (5-2) 


where the second equality follows from the fact that Ti = T A ri, and ri is a random 
variable independent of X, with density function /3e“'^^l{j>o}. 

(ii) Next, notice that for any t and bounded continuous function (po, one has by 
integration by parts that 


dxE[Mx + bt + Wt)] 


E 


(j)o{x + bt + Wt) 


m- 

t. 


(5.3) 


Notice also that ARti = Wt^, ATi = Ti and Xt^ := xq + bTi + Rti- It follows by 
Lemma A.3 that 


dxu{0,xo) = E 




ARti 

ATi 


(^fl'(AT)l{T<Ti} + ^—^93;u(ri, Ari)l{Ti<T}) • (5.4) 


Changing the initial condition (0,xo) to {Ti,Xt^), one obtains that, whenever Ti < T, 


dxu(Ti, Xti) = E 


' ffATa AVErg / 


Ti,Xt, 


Plugging the above expression of dxu{Ti, Xti) into the r.h.s. of (5.2), and using the 
fact that T < r 2 is equivalent to Nt- < 1, and P[{A^ 7 ’_ < 1} \ {Nt < 1}] = 0, it 
follows that u(0 ,Xq) = E['i/;i]. 
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(iii) Next, changing the initial condition in (5.4) from (0, xq) to (T 2 , when T 2 < T, 
and then plugging the corresponding expression of dxu(T 2 , XT 2 ) into ijji, it follows that 
tt(0,xo) = lE['i/; 2 ]- Repeating the procedure, we have for all n > 0, 

E[g{XT)] = u{0,xo) = E['i/’n]- 


Finally sending n —>■ 00 , and using Lemma A.2 together with the dominated conver¬ 
gence theorem, it follows that E[( 7 (A'r)] = E[lim„^oo V'n] = E['i/;]. □ 


Remark 5.2. We can also interpret formally the representation ip in (5.1) as the 
expansion of the diffusion process X around a Brownian motion. Let 6 = 0 and 
p.{t,x) = go for some constant go S M, so that Xt = Wt and Xt = got + Wt- Using 
the fact that ¥{Nt = k) = e~^'^ , V/c > 0, it follows formally that 


E[g{XT)] =E 




Lfc=0 


k\ 


= E 




where the second equality follows by the fact that 


(5.5) 


r[NT = k]= 


k\ 


and E 


-,{k) 




i + l 


2 = 0 


= E 


K 


TT^ 


i=0 


AWi+1 

AT,- 


i+i 


In particular, the r.h.s. of (5.5) is exactly defined by (5.1) in this case. 

To conclude this part, we notice that ip in (5.1) is integrable, but has an infinite 
variance in general. In the next subsection, we exploit the arbitrariness of the constant 
6 which is involved in the definition of ip. More precisely, we shall choose different 
constants 6 at each time T^, in an adaptive way. This will lead to the estimator pj in 
Theorems 2.2 and 3.3. 


5.2 A regime switching diffusion representation 

For d > 1, T > 0, let {g, a) : [0, T] x —)• x be bounded continuous functions 

satisfying 

\g{t,x) - g{t,y)\ + \a{t,x) - a{t,y)\ < L\x - yp (t,x,y) € [0,T] x M"* x (5.6) 

for some constant L > 0. We start by considering a linear parabolic PDF 

dtu + g ■ Du + a : D^u = 0, on [0,T) x M'^, (5.7) 

with terminal condition u{T,x) = g{x), where a(-) := ^crcr"'"(•), A : B := Tt{AB~^) 
for any two d x d dimensional matrices A,B ^ M'^, and D, denote the gradient 
and Hessian operators with respect to the space variable x. Next, let us consider the 
diffusion process (A's’^°)<je[o,T] defined as unique strong solution of the SDE 

Aq = xo, and dXg = g(^s,Xs) ds + a[s,Xs') dWg, sG [0,T]. (5.8) 

When PDF (5.7) admits a classical solution in C'j[’^([0, T] x W^), i.e. the collection of all 
functions (p{t, x) such that (p, dtf, Dcp, D'^cp and D^cp all exit and bounded continuous, 
it follows by Feynmann-Kac formula that Vq := E[gf(A.^’^°)] = u(0, xq). 
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Remark 5 . 3 . For technical reason, we will assume that u € Cl'^{[0,T] x rather 
than in T] X M'^). But by approximating the coefficient pL,a with smooth func¬ 

tions, one can relax this regularity condition in more concrete context. 

Recall that for /3 > 0, (Ti)j>o is a sequence of i.i.d. £’(/3)-exponential random 
variables, which is independent of the Brownian motion W. We define 

k 

Tk := A: > 0, and Nt := max |/c : T^. < t}. 

i=l 

Then {Nt)o<t<T is a Poisson process with intensity (3 and arrival times {Tk)k>o^ and 
To = 0. We also introduce, for all k > 0 , AWj^ := bP(Tfe_i+t)ATfe - It is clear 

that the sequence of processes (AW.^)fc>o are mutually independent. 

Let {ft, d) : {s, y, t, x) € [0, T] x x [0, T] x x be uniformly bounded, 

and continuous in t, Lipschitz in x, we define X by 

Xo := xq and dXt = il{Bt,t, Xt)dt + Xt)dWt, (5.9) 

with Qt := (T/Vt, ). In other words, the process X is defined recursively by, 

Xq = xq and for all A: > 0, 

XTk+i=^Tk + / il[Tk,XT^,s,Xs)ds + / d[Tk,XT,^,s,Xs)dWs. 

Jn Jn 

Example 5 . 4 . (i) Let {fi,d){s,y,t,x) = {y,a){s,y), then X is defined as a Euler 
scheme as in (3.2), i.e. Xq = xq, and 

= Xt^+ h{Tk,XTjATk+i + a{Tk,XTjAWT^,+.^. 


(ii) When /!{■) = 0 and d{s, y, t, x) = a(s, y) + dxiy{s, y){x — y), then SDE (5.9) turns 
to be a linear SDE, whose solution is given explicitly in (3.8). 


We first formulate an assumption on the existence of automatic differentiation 
weights associated to SDE (5.9). Let 6 G [0,T) x and (t,x) G [0, T] x the 
process (-^s’^’^)se[t,r] is defined by SDE 

Xl'^'^ :=x, dXi’^’^ = fi{e,s,Xl’^’^)ds + d{e,s,Xi’^’^)dWs, (5.10) 

Assumption 5.5. There is a pair of measurable functions (W0(‘))I^('))7 called auto¬ 
matic differentiation weights, taking values mM'^xM'^, such that, for all 9 ^ [0, r)xM'^, 
(t, x) G [0, T) X s > t, one has (LVg(A, x,s — t, {Wr — Wt)r^[t,s\)^ * = 1)2) are both 
integrahle. Moreover, for all bounded continuous function —)• M, 




* = 1 , 2 , 


where D,D^ denote the gradient and Hessian operators with respect to the variable x. 
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Let a(-) := ^aa~^{■) and a(-) := ^a'a~^{■), we denote 

00 = {t,x) and then 0 ^ = {TI,Xt^.), for all k > 0. 


and then for k > 0, 

Afk := {fi,a){n,XT,)-{fi,a){&k-i,Tk,XT,) G 


and for A: > 0, 

Wfc := (wi^,i^J(rfc,XT,,rfc+i,Aiy.'^+i) g 

with the weight functions (W’ 0 (-), (•)) given in Assumption 5.5. We then define 

Nt 

^ := (<7(Xt)-5(^t^^)1{7V^>o}) (5.11) 

k=l 

where {p, P) • {q,Q) := p ■ q + P : Q for all p,q € P,Q & M*^. Here we use the 
convention H^^j^ = 1. Finally, for all n > 1, we also introduce 


NtAu 


ipn = e- 


— ^/^^n + 1 


n (r'A/fc.Wfc) UXt) - g{XT^^)l 


k=l 


{Nt> 0 } ) ^{NT<n} 


+ /3 ^(A/n+i« {Du, D^u)(Tn+l,XT^_^_^)')l{NT>n}{^A2) 


Assumption 5.6. (i) The sequence {'ipn)n>o is uniformly integrable. 

(ii) Let {ei}i=i^... denote the canonical basis o/M'^. There is some £q > 0, such that 
for all {t, x) G [0, T) x and 6 G [0, T) x n > 0 and i = 1, - ■ ■ ,d, all the following 
random vectors is integrable: 


W^{t,X,Ti A (T - t), {Wr - Wt)r(z[t,{t+ri)AT]), 


sup 

ee( 0 ,£o] 


Wl {t, x + £ei,TiA{T-t), {W. - Wt)) - Wl {t, X, 


TiA{T-t),{W.-Wt)) 


and 

Afn+i • (Hu,D^n) (r„+i, AT^+JWn. 


Theorem 5.7. Suppose that the PDE (5.7) has a classical solution u G {[0,T] x 
suppose in addition that Assumptions 5.5 and 5.6 hold true. Then if is integrable 
and tt(0, xq) = E [')/)]. 


Remark 5.8. (i) The condition that u G C'^’^([0,T] x may be relaxed in the 
concrete applications of Theorem 5.7. This will be indeed performed in Section 3.3 
by exploiting the integrability of the automatic differentiation weights (LVg,VV|) of 
Assumption 5.5. 

(ii) By definition, the automatic differentiation weight satisfies E[VVfc] = 0, then if in 
(5.11) has the same mean than the estimator 


g{XT) 


Nt 

n (A/fc.Wfc). 


k=l 
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However, in practice, the weight function Wfc is typically of inhnity variance, or even 
not integrable, in general. Indeed, as we will see in the following, Wfc is generally 
of order where conditioning on Nt = n, (Ti, • • • ,TjVj,) follows the 

law of statistic order of uniform distribution on [0, T]. Then by direct computation, 
one knows E[l/ATjVj,_|_i] = oo. In the dehnition of 'i/’ in (5.11), the additional term 
l{jVj.>o} can be seen as a control variate so as to guarantee the integrability 

of 'll). 

(iii) As a consequence of the integrability problems raised in (ii). Assumption 5.6 is 
in fact implicitly a restriction on the choice of the coefficients jl and a, and we cannot 
expect a representation for u{t, x) with arbitrary fi and a, see Section 5.3 below. 

Proof of Theorem 5.7. (i) Recall that u € T] x M'^) is a classical solution 

of PDE (5.7). Denote {fi 0 ,ao){-) = one can rewrite (5.7) in the following 

equivalent way: 

—dtu — fie ■ Du — ag : D^u — ((// — fig) ■ Du + (a — ag) : D^u) = 0 . (5.13) 


Using Feynmann-Kac formula, it follows that 

u{0,xo) + [ {{g-flg)-Du + {a-ag):D‘^u){s,X^'^°’‘^)ds\, (5.14) 

L Jo 

where is dehned by (5.10), which coincides with X in (5.9) on [0, Ti] whenever 

0 = (0,xo). 

Recall that Ti = ri A T, where ri is a random variable of density 
independent of the Brownian motion W. Fixing 9 = (0,xo), it follows that 

u(0, xo) = E(At) l{Arj,=o} +/3"^A/i • (Du, D^u)(Ti, Ati) l{Ar,j,>o} 


(ii) Let us now go back to the expression (5.14), and derive an expression for the 
derivatives Du{0,xo) and D‘^u{0,xo). First, for Du(0, xq), we use the integrability 
condition in Assumption 5.6 with Lemma A.3, and also the fact the Du{-) is continuous, 
it follows that 


Du{0, xq) 


E 


g{X^^°’^)W^g{xo,T) 


+ 


/ {{g - fig) ■ Du + {a 

Jo 


dg) : D‘^u){s,X^’^°'^)Wg{xo,s)ds , 


where we simplify the notation (0, xq, s, (Wr — ^t)re[o,s]) ^^(^^Oj'S)- Then by 

the independence of ri to the Brownian motion W, and setting 6 = (0,xo), it follows 
that 


Du(0,xo) = E -00 >V(o,^g)(0,xo,ri, AIT.^) 


(5.15) 


Next, for D^n(0, xq), we use again Lemma A. 3 together with the integrability 
condition in Assumption 5.6 and Lipschitz property of x i-A ((^ — /i) • Du + (a — dg) : 
D'^u){s, x), and the continuity of D^u(-) that 


dVo,xo) = d 2 ^E[ 5 (A°’"°’^)]+/ D^^E 

40 


{{g-fig)-Du+{a-dg) : D^u) (s, 


ds. 
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Setting 0 = (0, xo) and using Assumption 5.5, it leads to 


D‘^u{0,xo) = E ^0 >Vfo,2,Q)(0,xo,ri, 

Recall that E[>V^^g xq, T, AVh.)] = 0, we then obtain 

D^u{0,xo) = E {ipo - e^^g{xo)l{NT=o}) ^o,xo)(0> ®o, AIR.^) 


(5.16) 


(iii) Changing the initial condition (0,xo) in (5.15) and (5.16) by (Ti, Ai) (remember 
^jJQ dependent also on the initial condition (0,xo)), then plugging the expression of 
Zl^u(Ti, Ai) and D'^u{Ti, Xi) into the definition of ipo in (5.12), it follows by identifying 
the term that 


u(0,xo) = E[V^i]. 

(iv) Repeating the arguments by replacing the initial condition (0, xq) by (T„+i, A„+i) 
in (5.15) and (5.16) and then plugging the corresponding expression into the dehnition 
of V’nj etc., we obtain that tt(0,xo) = E['i/;„] for all n > 0. Then letting n —> oo, we 
obtain 


u(0,xo) = lim E[^n] = E[ lim ijjn] = 

which concludes the proof. □ 

5.3 Proof of the representation results in Theorems 2.2, 
2.5 and 3.3. 

Using the results in Theorem 5.7, we can easily complete the proof of the representation 
results in Theorems 2.2, 2.5 and 3.3. 

Proof of Theorems 2.2 (ii) and 3.3 (ii). (i) In the context of Theorems 2.2 and 3.3, 
the increment — Xt^, conditional on (T^, At^^.), is Gaussian. And the estimator 

Tp corresponds to the estimator in Theorem 5.7 with automatic differentiation weights 
function 

Wg{-,St,Sw) := and W^{-,St,Sw) ■= ^(5.17) 

In particular, it is clear that Assumption 5.5 holds true with the above choice of 
automatic differentiation weight functions in (5.17). 

(ii) Besides, the uniform integrability conditions and integrability conditions in As¬ 
sumption 5.6 can be easily obtained following the lines in the hrst part of the proof of 
Theorems 2.2 and 2.2, using Lemma A.2. 

(iii) Now, suppose in addition that g, a and g are bounded smooth functions with 
bounded continuous derivatives, so that u £ C'^’^([0, T] x M'^). It follows by Theorem 
5.7 that Vo = E[V']. 
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(iv) Finally, when fi(-) and cr(-) satisfy the Lipschitz condition (3.1) and g is Lipschitz, 
we can find a sequence of bounded smooth functions g^{-)) which converges 

locally uniformly to g[-)) as e —>■ 0. Let X'^ be the solution of 

dXt = ge{t,Xf)dt + a,it,X^)dWt. 

Then by the stability of SDEs together with dominated convergence theorem, it follows 
that 


:= E[g,{X^)] Vo :=E[g{XT)], as e ^ 0. 

Moreover, by Lemma A.2 together with dominated convergence theorem, it is easy to 
prove that E['i/’^] E[^] as e ^ 0, where denotes the estimator of the algorithm 

(3.3) associated to the coefficient ge)- We then conclude the proof. □ 

Proof of Theorem 2.5 (ii). For the path-dependent case, it is enough to use the 
same arguments as in Theorem 2.2, together with the PDF system (2.12) in Remark 
2.3. □ 


5.4 Proof of Theorem 3.5 

To introduce the algorithm in the context of Theorem 5.7, we propose to choose 
/!(•) = 0 and d{s,y,t,x) = a{s,y) + d^a{s,y){x - y). 


Before providing the proof of Theorem 3.5, we first give a lemma which justihes 

_2 

our choice of the automatic differentiation weight function in (3.11), as well as 
some related estimations. Let ci,C 2 ,x G M be constants such that ci -|- C 2 X ^ 0, we 
denote by X^'^ solution of the SDE 

Ao = x, dXt = [ci + C2Xt)dWt, (5.18) 


whose solution is given explicitly by 


_ 

Aj — < 


“ § + (§ + , if C2 / 0, 


X ciWt, 

Consider also its antithetic variable Xf defined by 


= 


X - ciWt-, 


if C2 = 0. 


“ § + (§ + ^) , if C2 7^ 0, 


if C2 = 0. 


(5.19) 


Lemma 5.9. Let x G M, (ci,C 2 ) G he two constants such that ci + C 2 X 7 ^ 0, 
(/> ; M —)• M a bounded continuous function. 

(i) Then for all t £ {0,T], 


CE[<(<(X?’^)] 



1 

(Cl -I- C2 x)2 


Wt , 

-C2- + 


f^ )v 


(5.20) 
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(ii) Suppose in addition that Then there is some constant C independent 

of {t, x) such that, for all {t, x) G [0, T] x 


E 


xn - 




2l 


+ E 




0,x 


< C(ci + C 2 xf 


t2 ) \ 


Proof, (i) First, when C 2 = 0, it is clear that result is correct (see e.g. Lemma 2.1 of 
Fahim, Touzi and Warin [13]). Next, when C 2 / 0, denote v{x) := then 

with the expression of Xf in (5.19), it follows that 

v(x) = [ (!>{ — — + {— + a;')e“'^ 2 V 2 +c 2 v^y'\ 

Jk. ^ C2 Vc2 / / V27r 

Suppose that </>(•)£ then using integration by parts, it follows that 


v'{x) 


I 

Jr 

j 

Jr 


- — + [ — + x]e 
C2 


- — + 
C2 


Cl 

C2 

Cl 


-c|t/2+C2\/iy'\ -c^t/2+C2Viy. 




e 




C2 


E 


cPixX 


1 


+ x]e 

W, 


3-c|l/2+C2\/tj/ 




Cl + C2X i/t 




Cl + C2X t 


Similarly, still using integration by parts, and by direct computation, we obtain 


v"(x) 


E 




(ci + C 2 X)^ 


Wt , 
-C2- + 


WXt ^i 

t^ ). ■ 


When (/>(•) is only a bounded continuous function, one can approximate by a 
sequence of smooth function (fsi') which converges to 4>{-) uniformly, and (f's and <f>f 
are bounded continuous. We then obtain 


v^(x) := E[(j)s{xX] v{x). 


Moreover, the limit lim^^o c(,(a:), lim£^oC£(x) exist, thus v"{x) also exists and 


v"(x) 


lim vf(x) 
s^o ’ 



1 

(Ci + C2 x)2 


Wt , 

-C2- + 


wj-u i 

t^ ). ■ 


(ii) When C 2 = 0, the estimation in (ii) of the statement is clear true since (j)' and cj)" 
are uniformly bounded. 

When C 2 / 0, denote j^'jo := sup^. \(j)'{x)\, we obtain, by direct computation, that 


E 





< 


|</.'|oE (X 


J,ir 



2Xl\ 

. 


= |(/>'|oE (ci+C2x)' 


^-clm+c^Wt _ l.2W‘l{c2Wt - clt/2f^ 
C 2 Wt - clt/2 ) t 2 


which is clearly uniformly bounded by C(ci + C 2 x)^ for some constant C independent 
of {t, x) G [0, T] X 
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Next, denote i{y) := {x + (e '^ 2 */ 2 +c 2 !/ _ define ip{y) := (pi^ + ^{y))- 

Then 

^''{y) = </."(x + ^(y))(c2 + cix)2e-^2i+2c23/ + + cix)c2e-^=*/2+'^=2^.(5.21) 


It follows by the definition of as well as its derivative, together with direct compu¬ 
tation, that 


E 

= E 
< E 


{piWt) + pi-Wt) - 2c/p(0)) 


O.x 


2/Ilf -tN2 


t2 


-h E 


„/W2_i,2 


t2 


t2 


sup p"{z) 
\z\<\Wt\ 


+ E 


X 




which is also uniformly bounded by C{ci -|- C 2 x)^ for some constant C > 0 , □ 

Proof of Theorem 3.5. (i) Let us first prove that E[^^] < oo for 'tp dehned by 
(3.12). First, we notice that = W| for all /c = I,-- - ,Nt — 1, 5 S (72(1^), and 
with the choice of c\ and C 2 in (3.7), one has c\ + = a{Tk, XTk)j which is 

uniformly bounded. By considering the conditional expectation over (X^jv^, ATat^+i) 
using items (ii) of Lemma 5.9, we have E[|'0|^] is bounded by 


CE 


-2Nt TT / Xt^P) — ak 
2a{Tk,XTfP) 


/3- 


n{ 

k=2 


„ _ ^ AWT.^Wl-ATk 

— dx(y{Tk,XTfP) -I- 


An 




)n 


for some constant C. Further, by denoting AXt^. '■= Xxf, — Xt^_^, one has 

|a(7fc, Atj.) — afc| < (^\a\o + \dxCr{Tk-i, XT^_-^)AXT,,\/2'j 

(\dta\oAn + \d‘^xx^l{AXT,f), 

where |cr|o := sup^ x ^)l- Notice that a > £ > 0, a and dxcr are uniformly bounded, 
then to prove that ip is of hnite variance, it is enough to prove that, for some (7 > 0 
large enough, the expectation of 


Nt 

n [C(C +\dxa{n-uXT,_,)AXT, 


k=2 


c + 


AXl 

An 


\c\AWt, 


+ 


An 


+ 1 


2i 


(5.22) 


is hnite. Similarly to the computation in item (ii) of Lemma 5.9, we have 


AAt, = Xt,-Xt,_, 

exp ( - dxa{n-i,XT,_,fAn/2 + dxa{n-i, Xt,_,)AWt,) - 1 


= a{n-l,XT^:_l)- 


dx<^{Tk-i.,Xk-i) 
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Notice again that cr(-) and dxcr^-) are uniformly bounded, it follows that 


e{ [(c + \dxa{Tk-uXT,_,)AXT,\) (C + (c| Abhr, 


+ 


AW2 

_£fc 

ATfc 


+ 1 


1 2 


XT,_„Tk-i,ATk] < C", 


for some constant C' > 0 independent of Xj'f ._^, Tfc-i) Then the variance of (5.22) 
is bounded by CE[(C')'^^] < oo and hence ifj in (3.12) is of hnite variance. 

(ii) Let us now consider the estimator 'tp. By the same computation, we obtain that 

E[ ^ I A't, ATi, • • • , ATjVj,+i] < —= ^ , for some C > 0, 

V ATvy+i 

where the r.h.s. is integrable but of inhnite variance (see Lemma A.2). Similarly, it is 
easy to check the uniform integrability condition in Assumption 5.6 for ?/) in (3.10). 

(iii) Finally, using Lemma 5.9 (i), it follows that Assumption 5.5 holds true. Moreover, 

with the regularity condition on x) and g in Assumption 3.4, we know u € C'^’^(M). 
We then conclude the proof of u(0, xq) = E[V’] = E[?/)] by Theorem 5.7. □ 


A Appendix 

We first provide an estimation on the order statistics of uniform distribution on [0,1], 
which induces an estimation on a functional of the arrival times {Tk)k>o of the Poisson 
process. We next provide a technical result on the automatic differentiation function 
related to a SDE. 


Lemma A.l. Let p € (0,1), {Uk)k=i,- ,m a sequence of i.i.d. random variable of 
uniform distribution on [0,1], and (f7(i) < t7(2) < ••• < t7(m) ^he associated order 
statistics. Then 


G 


m,p 


E 


(— --- 


^(m) ^{m—1) 


< ml 


(i-pY 


Proof. First, we notice that for any x G (0,1), 



(A.l) 


Then, since the density of the order statistics (77(1), • • • , t7(^)) is provided by 


f(ui,--- ,Um) ■ ml 


it follows by direct computation that 


G 


m 


= ml 



1 


1 


u^-l ^UiU2- Ui 



• • • dUm 


< ml 


1 

(1 -p)"^ 


where the last inequality follows from (A.l). 


□ 
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Let N = {Ns)s>o be a Poisson process with arrival times {Tk)k>o, denote AT^+i := 
Tfc+i — Tfc. Let 0 = to < ti <■■■< tn = T < oo be a discrete time grid, we define 
further 

Tfc := mm(Tk,ti), whenever Tj^-i G [ti-i,ti) for some i = 1, • • • ,n, 


and 


ATfc := fk - Tfc_i for every k = 2,--- ,Nt + 1. 


Lemma A.2. Let p G (0,1), then for every constant C > 0, one has 


n 


C 1 

(ATfc+OP - 


< oo. 


Proof, (i) Notice that we can always add points into the time grid 0 = to < < 

' ' ' tn = T, which makes ATj+i smaller. Therefore, one can suppose without loss of 

generality that tk — tfc-i < (2/3(7)^“^ for every /c = 1, • • • , n. 

(ii) For every A; = 1, • • • , n, we denote := ff{i : Ti G [tk-i,tk)}, and := 

ffc, with ki := + * for i = I,-- - , + 1, and ATi^ := By 

the memoryless property of the exponential distribution, it is clear that (AT-^, i = 
2, • • • + l), • • • , (AT”, i = 2, - ■ ■ A” + l) are mutually independent. Moreover, we 

have 


Nt ^ 

TT 

(AT,+i)P 


n / N'^ 

n n 

k=i \i=i 



(A.2) 


Next, the law of [Th, i = I,-- - ,A^) conditioning on A^ = m is the law of order 
statistics of uniform distribution on [tk-i,tk]- Then it follows by Lemma A.l that for 
every A; = 1, • • • , n, 


E 


C 


AT'' 


g Pij'k ^k — l) 


iP{tk -tk-i)y 


m=0 


-m!2" 


ml 


C 


{tk-tk-i)P 


m=0 


where the last inequality follows by the fact tk — Afc-i < (2/3(7)^ We then conclude 
the proof by (A.2). □ 


Let be the solution of SDE 


XS = x, dXf = fi{t,Xf)dt + a{t,Xf)dWt, 

where {p, a) : [0, T] x —)• x M'^ is continuous and in addition Lipschitz continuous 

in X. 
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Lemma A. 3. Suppose that for all t £ [0, T] and bounded continuous function cj) : 
[0,r] X M, the derivatives [dxi^[4>{t,Xf)^,d‘^.^^M[(l){t,Xf)]).^^-^ ^ exist; and 

there is some measurable W^-valued function W^[x,t, (W4)sg[o,t]) such that for 

dxM[cfit,Xf)] = E[(P{t,Xf)Wlix,t,{Ws)s^io,t])], i = l,---,d. 

Let F{dt) be some probability measure on [0, T]. 

(i) Suppose that for each x € i = 1, - ■ ■ ,d, 

rT 


Jo 


< oo, 


and the continuous function (j) : [0, T] x ^ M (which may be unbounded) satisfies 


E 


\f)it,Xf)W}{x,t,{Ws)s^[o,t])\\F{dt) < oo. 


Then 


dx, rE[<fit,Xf)]F{dt) = riKU{t,Xf)Whx,t,iWs)sem)\d^idt). (A.3) 
Jo Jo 


to ' " Jo 

(ii) Suppose in addition that (j){t, x) is bounded continuous function and Lipschitz in 
x, and for each x £ i,j = l,---d, 
rT 


F{dt) < oo, 


and 

r sup 

J 0 £■£ [0,£o] 


-E 


w }(x + ecj , t, (W's) 56 [ 0 ,t] )-yVl{x,t, (TVs) 56 [ 0 ,t]) F{dt) < oo, (A.4) 


for some Eq > 0, where [ej)j=i^...^d denotes the canonical basis o/M“. Then 

l\[cP{t,Xf)]F{dt) = £ dl,M[c){t,xn]F{dt), (A.5) 

where, in particular, the partial derivative at the l.h.s. and the integration at the r.h.s. 
are well defined. 

Proof, (i) First, let us notice that {t,x) e-)- {pL,a){t,x) is Lipschitz in x, then by 
standard analysis (see e.g. Chapter 7.8 of [17]), there is some constant C independent 
of e > 0 and i = 1, ■ ■ ■ d, such that 

J^X+EBi _ J^X 2^ 


E 


< (7(1 + 6' 


Ct\ 


(A.6) 


(ii) Suppose that (f(t, x) is bounded continuous and Lipschitz in x. It follows that 

rT 


hm / -E\{(){t,Xf^^‘^^)-),{t,Xn)]F{dt) 


= E 


[ - ^{t,Xf))]F{dt) 

r <Pit,Xf)W}{x,t, iWs)se[o,t])F{dt) 

Jo 
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where the first equality follows by the Lipschitz property of x (p{t, x) and (A.6). We 
hence proved (A.3) when x >->■ x) is Lipschitz. 

(hi) When (j) is only continuous, it is enough to approximate it by a sequence {4>n)n>i 
which are all bounded, and Lipschitz in x. Then by the integrability of >V/(x, t, 
as well as that of 4>{t, Af )>V/(x, t, (Ws) 5 g[o,t]) under F{du}) x F{dt), it follows that (A. 3) 
holds true for continuous function cj). 

(iv) To prove (A.5), let us use (A.3) and obtain that 


lim — 
£—^0 S 


d^. r ¥.[ct>{t,X^+^^^)]F{dt)-d^^ r ¥.[<i>{t,Xf)\F{dt) 
Jo Jo 

r 

/ -E </>(L- 

Jo ^ L 

F 1 r 

/ lim -E (bit,. 

Jo ^ ^ 


= lim 
£—^0 


Xj+^^^)WXx + eei,t,-) - (j){t,Xf)WUx,t,-) F{dt) 


{x + eei,t,-) - ^{t,XJ)W {x,t,-) F(dt). (A.7) 


where the first equality follows by the Lipschitz property of x eT” 4>{t, x) and the 
estimation (A.6) together with (A.4), and in particular, the integrable in the last term 
of (A.7) is well defined, and hence the limit of the first term of (A.7) exists. □ 
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